L’objectif de ce TP est d’illustrer les notions abordées dans le chapitre de régression linéaire.
Les librairies R nécessaires pour ce TP :
library(ellipse)
library(leaps)
library(MASS)
library(corrplot)
library(glmnet)
library(coefplot)
library(ggplot2)
library(gridExtra)
library(ggfortify)
library(plotly)
library(reshape2)
La pollution de l’air constitue actuellement une des préoccupations majeures de santé publique. De nombreuses études épidémiologiques ont permis de mettre en évidence l’influence sur la santé de certains composés chimiques comme le dioxyde souffre (SO2), le dioxyde d’azote (NO2), l’ozone (O3) ou des particules en suspension. Des associations de surveillance de la qualité de l’air (Air Breizh en Bretagne depuis 1994) existent sur tout le territoire métropolitain et mesurent la concentration des polluants. Elles enregistrent également les conditions météorologiques comme la température, la nébulosité, le vent, les chutes de pluie en relation avec les services de Météo France… L’une des missions de ces associations est de construire des modèles de prévision de la concentration en ozone du lendemain à partir des données disponibles du jour : observations et prévisions de Météo France. Plus précisément, il s’agit d’anticiper l’occurrence ou non d’un dépassement légal du pic d’ozone (180 \(\mu\)gr/m3) le lendemain afin d’aider les services préfectoraux à prendre les décisions nécessaires de prévention : confinement des personnes à risque, limitation du trafic routier. Plus modestement, l’objectif de cette étude est de mettre en évidence l’influence de certains paramètres sur la concentration d’ozone (en \(\mu\)gr/m3) et différentes variables observées ou leur prévision. Les 112 données étudiées ont été recueillies à Rennes durant l’été 2001.
Les 13 variables observées sont :
Les données sont disponibles sur la page moodle du cours. Pour les importer, vous pouvez utiliser la commande suivante :
ozone = read.table("Ozone.txt")
Afin de vous familiariser avec les données, faites dans un premier
temps une analyse de statistique descriptive. Vous pouvez utiliser les
fonctions summary(), boxplot(),
pairs(), barplot(), corrplot(),
….
# A COMPLETER - Faire des stat descriptives des données
summary(ozone)
maxO3 T9 T12 T15
Min. : 42.00 Min. :11.30 Min. :14.00 Min. :14.90
1st Qu.: 70.75 1st Qu.:16.20 1st Qu.:18.60 1st Qu.:19.27
Median : 81.50 Median :17.80 Median :20.55 Median :22.05
Mean : 90.30 Mean :18.36 Mean :21.53 Mean :22.63
3rd Qu.:106.00 3rd Qu.:19.93 3rd Qu.:23.55 3rd Qu.:25.40
Ne9 Ne12 Ne15 Vx9
Min. :0.000 Min. :0.000 Min. :0.00 Min. :-7.8785
1st Qu.:3.000 1st Qu.:4.000 1st Qu.:3.00 1st Qu.:-3.2765
Median :6.000 Median :5.000 Median :5.00 Median :-0.8660
Mean :4.929 Mean :5.018 Mean :4.83 Mean :-1.2143
3rd Qu.:7.000 3rd Qu.:7.000 3rd Qu.:7.00 3rd Qu.: 0.6946
Vx12 Vx15 maxO3v vent
Min. :-7.878 Min. :-9.000 Min. : 42.00 Length:112
1st Qu.:-3.565 1st Qu.:-3.939 1st Qu.: 71.00 Class :character
Median :-1.879 Median :-1.550 Median : 82.50 Mode :character
Mean :-1.611 Mean :-1.691 Mean : 90.57
3rd Qu.: 0.000 3rd Qu.: 0.000 3rd Qu.:106.00
pluie
Length:112
Class :character
Mode :character
[ getOption("max.print") est atteint -- ligne 1 omise ]
Dans cette section, nous souhaitons expliquer la concentration d’ozone maximale de la journée (maxO3) en fonction de la température à 12h (T12).
plot() (ou geom_point() de
ggplot2).# A completer
# Scatter plot
plot(ozone$T12, ozone$maxO3,
xlab = "Temperature at 12:00 PM (T12)",
ylab = "Maximum Ozone Concentration (maxO3)",
main = "Scatter Plot of maxO3 vs. T12")
Question : Pensez-vous que l’ajustement d’un modèle de régression linéaire \(y_i = \theta_0+\theta_1 x_i +\varepsilon_i\) est justifié graphiquement?
Réponse : On voit une espèce de tendance linéaire : un modèle de régression linéaire est justifiée.
lm() et exploitez les résultats.# Régression linéaire simple
reg.simple = lm(maxO3 ~ T12, data = ozone)
# Afficher un résumé des résultats de la régression
summary(reg.simple)
Call:
lm(formula = maxO3 ~ T12, data = ozone)
Residuals:
Min 1Q Median 3Q Max
-38.079 -12.735 0.257 11.003 44.671
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -27.4196 9.0335 -3.035 0.003 **
T12 5.4687 0.4125 13.258 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 17.57 on 110 degrees of freedom
Multiple R-squared: 0.6151, Adjusted R-squared: 0.6116
F-statistic: 175.8 on 1 and 110 DF, p-value: < 2.2e-16
Interprétation :
Que contiennent les sorties suivantes :
reg.simple$coefficients
reg.simple$residuals
Réponse : coefficients -> vecteur des valeurs observées de \(\^{\theta}^{obs}\) residuals -> vecteur des résidus
ggplot(ozone, aes(T12, maxO3))+
geom_point() +
geom_smooth(method=lm, se=TRUE)+
xlab("T12")+ ylab("maxO3")
Ce graphique permet visuellement de vérifier l’ajustement des données au modèle de régression linéaire. Que remarquez-vous?
autoplot(reg.simple,which=c(1,2,4),label.size=2)
On prendra soin de comprendre les différents graphiques obtenus.
temp_var = predict(reg.simple,
interval="prediction")
new_df = cbind(ozone, temp_var)
ggplot(new_df, aes(T12, maxO3))+
geom_point() +
geom_line(aes(y=lwr), color = "red",
linetype = "dashed")+
geom_line(aes(y=upr), color = "red",
linetype = "dashed")+
geom_smooth(method=lm, se=TRUE)+
xlab("T12")+
ylab("maxO3")
confint(), construisez un
intervalle de confiance pour chaque paramètre séparément.# A COMPLETER
# Intervalles de confiance pour theta_0 (intercept) et theta_1 (pente)
IC = confint(reg.simple)
IC
Pour tenir compte de la dépendance entre \(\theta_0\) et \(\theta_1\), on peut aussi construire une région de confiance pour le vecteur \(\theta=(\theta_0,\theta_1)'\) grâce aux commandes suivantes. Comparez les résultats.
df1 = as.data.frame(
rbind(coefficients(reg.simple),
ellipse(reg.simple,level=0.95)))
colnames(df1) = c("intp", "slope")
ggplot(data=df1[-1,],aes(x=intp, y=slope))+
geom_path()+
annotate("rect",xmin=IC[1,1],xmax=IC[1,2],
ymin=IC[2,1],ymax=IC[2,2],fill="red",alpha=0.1)+
geom_point(data=df1[1,], aes(x=intp, y=slope),
pch=3)
Dans cette partie, nous souhaitons analyser la relation entre le maximum journalier de la concentration d’ozone (maxO3) et
On va donc utiliser le sous-jeu de données suivant
ozone1 = ozone[,1:11]
et mettre en place un modèle de régression linéaire multiple.
# Afficher un résumé statistique des variables
summary(ozone1)
# Afficher des boîtes à moustaches pour visualiser la distribution des variables numériques
boxplot(ozone1)
# Afficher une matrice de nuages de points pour explorer les relations entre les variables
pairs(ozone1)
\[ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + ... + \beta_p X_p + \varepsilon \]
lm() et commentez les résultats (on appelle
reg.mul la sortie de R dans la suite).reg.mul = lm(maxO3 ~ T9 + T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx12 + Vx15 + maxO3v, data = ozone1)
resume.mul = summary(reg.mul)
summary(reg.mul)
Call:
lm(formula = maxO3 ~ T9 + T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 +
Vx12 + Vx15 + maxO3v, data = ozone1)
Residuals:
Min 1Q Median 3Q Max
-53.566 -8.727 -0.403 7.599 39.458
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.24442 13.47190 0.909 0.3656
T9 -0.01901 1.12515 -0.017 0.9866
T12 2.22115 1.43294 1.550 0.1243
T15 0.55853 1.14464 0.488 0.6266
Ne9 -2.18909 0.93824 -2.333 0.0216 *
Ne12 -0.42102 1.36766 -0.308 0.7588
Ne15 0.18373 1.00279 0.183 0.8550
Vx9 0.94791 0.91228 1.039 0.3013
Vx12 0.03120 1.05523 0.030 0.9765
Vx15 0.41859 0.91568 0.457 0.6486
maxO3v 0.35198 0.06289 5.597 1.88e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 14.36 on 101 degrees of freedom
Multiple R-squared: 0.7638, Adjusted R-squared: 0.7405
F-statistic: 32.67 on 10 and 101 DF, p-value: < 2.2e-16
autoplot(reg.mul,which=c(1,2,4),label.size=2)
Dans le summary(reg.mul), un test est fait sur chaque
coefficient. Il revient à tester que la variable n’apporte pas
d’information supplémentaire sachant que toutes les autres variables
sont dans le modèle. Il est donc préférable d’utiliser des procédures de
choix de modèles pour sélectionner les variables significatives. On va
ici comparer la sélection de variable obtenue par différents critères:
BIC, AIC, \(R^2\) ajusté, Cp de
Mallows. Pour cela, vous pouvez utiliser la fonction
regsubsets() de la librairie leaps et la fonction
stepAIC() de la librairie MASS. Commentez les
résultats obtenus avec les différents critères, vous pourrez vous aider
des commandes suivantes :
# Utiliser regsubsets pour sélectionner les variables avec différents critères
choix = regsubsets(maxO3 ~ T9 + T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx12 + Vx15 + maxO3v, data = ozone1)
plot(choix, scale = 'Cp')
plot(choix, scale = 'adjr2')
plot(choix, scale = 'bic')
# Utiliser stepAIC pour sélectionner le modèle basé sur l'AIC
stepAIC(reg.mul)
Start: AIC=607.26
maxO3 ~ T9 + T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx12 + Vx15 +
maxO3v
Df Sum of Sq RSS AIC
- T9 1 0.1 20827 605.26
- Vx12 1 0.2 20827 605.26
- Ne15 1 6.9 20834 605.30
- Ne12 1 19.5 20847 605.36
- Vx15 1 43.1 20870 605.49
- T15 1 49.1 20876 605.52
- Vx9 1 222.6 21050 606.45
<none> 20827 607.26
- T12 1 495.5 21323 607.89
- Ne9 1 1122.6 21950 611.14
- maxO3v 1 6459.6 27287 635.51
Step: AIC=605.26
maxO3 ~ T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx12 + Vx15 + maxO3v
Df Sum of Sq RSS AIC
- Vx12 1 0.1 20827 603.26
- Ne15 1 6.9 20834 603.30
- Ne12 1 22.1 20849 603.38
- Vx15 1 43.5 20871 603.49
- T15 1 49.4 20877 603.52
- Vx9 1 264.0 21091 604.67
<none> 20827 605.26
- T12 1 641.4 21469 606.66
- Ne9 1 1183.6 22011 609.45
- maxO3v 1 7112.2 27940 636.16
Step: AIC=603.26
maxO3 ~ T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx15 + maxO3v
Df Sum of Sq RSS AIC
- Ne15 1 6.9 20834 601.30
- Ne12 1 23.1 20850 601.38
- T15 1 50.0 20877 601.53
- Vx15 1 79.6 20907 601.69
- Vx9 1 320.1 21148 602.97
<none> 20827 603.26
- T12 1 647.3 21475 604.69
- Ne9 1 1185.1 22013 607.46
- maxO3v 1 7112.6 27940 634.16
Step: AIC=601.3
maxO3 ~ T12 + T15 + Ne9 + Ne12 + Vx9 + Vx15 + maxO3v
Df Sum of Sq RSS AIC
- Ne12 1 16.2 20851 599.38
- T15 1 45.2 20880 599.54
- Vx15 1 75.1 20910 599.70
- Vx9 1 325.2 21160 601.03
<none> 20834 601.30
- T12 1 971.8 21806 604.40
- Ne9 1 1215.8 22050 605.65
- maxO3v 1 7106.3 27941 632.17
Step: AIC=599.38
maxO3 ~ T12 + T15 + Ne9 + Vx9 + Vx15 + maxO3v
Df Sum of Sq RSS AIC
- T15 1 45.6 20896 597.63
- Vx15 1 77.6 20928 597.80
- Vx9 1 337.6 21188 599.18
<none> 20851 599.38
- T12 1 1034.0 21885 602.80
- Ne9 1 2155.5 23006 608.40
- maxO3v 1 7130.2 27981 630.33
Step: AIC=597.63
maxO3 ~ T12 + Ne9 + Vx9 + Vx15 + maxO3v
Df Sum of Sq RSS AIC
- Vx15 1 74.1 20970 596.02
- Vx9 1 366.5 21263 597.58
<none> 20896 597.63
- Ne9 1 2241.1 23137 607.04
- T12 1 6712.6 27609 626.83
- maxO3v 1 7381.5 28278 629.51
Step: AIC=596.02
maxO3 ~ T12 + Ne9 + Vx9 + maxO3v
Df Sum of Sq RSS AIC
<none> 20970 596.02
- Vx9 1 903.4 21874 598.75
- Ne9 1 2714.8 23685 607.66
- T12 1 6650.4 27621 624.88
- maxO3v 1 7363.5 28334 627.73
Call:
lm(formula = maxO3 ~ T12 + Ne9 + Vx9 + maxO3v, data = ozone1)
Coefficients:
(Intercept) T12 Ne9 Vx9 maxO3v
12.6313 2.7641 -2.5154 1.2929 0.3548
Avec le critère BIC, nous retenons 4 variables : T12, Ne9, Vx9 et
maxO3v.
A l’aide des commandes suivantes, testez le sous-modèle avec les 4
variables retenues par BIC contre le modèle complet. Commentez.
reg.fin=lm(maxO3~T12+Ne9+Vx9+maxO3v,data=ozone1)
summary(reg.fin)
anova(reg.fin,reg.mul)
On commence par centrer et réduire les données avant de mettre en place et comparer des méthodes de régression régularisées.
tildeY=scale(ozone1[,1],center=T,scale=T)
tildeX=scale(ozone1[,-1],center=T,scale=T)
glmnet(), ajustez une
régression ridge en faisant varier \(\lambda\) sur une grille. On stockera le
résultat dans la variable fitridge. Explorez le contenu de
fitridge.lambda_seq<-10^(seq(-4,4,0.01))
# Ajuster la régression ridge avec glmnet
fitridge <- glmnet(x = tildeX, y = tildeY, alpha = 0, lambda = lambda_seq, family = c("gaussian"), intercept = F)
summary(fitridge)
Length Class Mode
a0 801 -none- numeric
beta 8010 dgCMatrix S4
df 801 -none- numeric
dim 2 -none- numeric
lambda 801 -none- numeric
dev.ratio 801 -none- numeric
nulldev 1 -none- numeric
npasses 1 -none- numeric
jerr 1 -none- numeric
offset 1 -none- logical
call 7 -none- call
nobs 1 -none- numeric
df=data.frame(lambda = rep(fitridge$lambda,ncol(tildeX)), theta=as.vector(t(fitridge$beta)),variable=rep(colnames(tildeX),each=length(fitridge$lambda)))
g1 = ggplot(df,aes(x=lambda,y=theta,col=variable))+
geom_line()+
theme(legend.position="bottom")+
scale_x_log10()
ggplotly(g1)
cv.glmnet() mettez en place une
validation croisée pour sélectionner le “meilleur” \(\lambda\) par MSE.# Réaliser une validation croisée pour sélectionner le meilleur lambda
ridge_cv <- cv.glmnet(x = tildeX, y = tildeY, alpha = 0, lambda = lambda_seq, nfolds = 10, type.measure = c("mse"), intercept = F)
# Trouver la meilleure valeur de lambda en minimisant le MSE
best_lambda <- ridge_cv$lambda.min
df2=data.frame(lambda=ridge_cv$lambda,MSE=ridge_cv$cvm,cvup=ridge_cv$cvup,cvlo=ridge_cv$cvlo)
gmse<-ggplot(df2)+
geom_line(aes(x=lambda,y=MSE))+
geom_vline(xintercept = ridge_cv$lambda.min,col="red",linetype="dotted")+
geom_line(aes(x=lambda,y=cvup),col="blue",linetype="dotted")+
geom_line(aes(x=lambda,y=cvlo),col="blue",linetype="dotted")+
#xlim(c(0,ridge_cv$lambda.min+0.5))+
scale_x_log10()
ggplotly(gmse)
La valeur de \(\lambda\) sélectionnée vaut
g1=g1 +
geom_vline(xintercept = best_lambda,linetype="dotted", color = "red")+
scale_x_log10()
glmnet(), ajustez une
régression Lasso en faisant varier \(\lambda\) sur une grille. On stockera le
résultat dans la variable fitlasso. Explorez le contenu de
fitlasso.# Ajuster la régression Lasso avec glmnet
fitlasso <- glmnet(x = tildeX, y = tildeY, alpha = , lambda = lambda_seq, family = c("gaussian"), intercept = F)
# Afficher un résumé des résultats
summary(fitlasso)
Length Class Mode
a0 801 -none- numeric
beta 8010 dgCMatrix S4
df 801 -none- numeric
dim 2 -none- numeric
lambda 801 -none- numeric
dev.ratio 801 -none- numeric
nulldev 1 -none- numeric
npasses 1 -none- numeric
jerr 1 -none- numeric
offset 1 -none- logical
call 6 -none- call
nobs 1 -none- numeric
df=data.frame(lambda = rep(fitlasso$lambda,ncol(tildeX)), theta=as.vector(t(fitlasso$beta)),variable=rep(colnames(tildeX),each=length(fitlasso$lambda)))
g3 = ggplot(df,aes(x=lambda,y=theta,col=variable))+
geom_line()+
theme(legend.position="bottom")+
scale_x_log10()
ggplotly(g3)
cv.glmnet() mettez en place une
validation croisée pour sélectionner le “meilleur” \(\lambda\) par MSE. En pratique, il est
préconisé d’utilisé lambda.1se (la plus grande valeur de
\(\lambda\) telle que l’erreur standard
se situe à moins de 1 de celle du minimum).lasso_cv <- cv.glmnet(...) # A COMPLETER
best_lambda <-lasso_cv$lambda.min
lambda1se <- lasso_cv$lambda.1se
La valeur de \(\lambda\) sélectionnée vaut ….
g3=g3 +
geom_vline(xintercept = best_lambda,linetype="dotted", color = "red")+
geom_vline(xintercept = lambda1se,linetype="dotted", color = "blue")+
scale_x_log10()
g3
Vous pouvez utiliser la fonction extract.coef() de la
librairie coefplot
# A COMPLETER
fitEN <- glmnet(...)
df=data.frame(lambda = rep(fitEN$lambda,ncol(tildeX)), theta=as.vector(t(fitEN$beta)),variable=rep(c(colnames(tildeX)),each=length(fitEN$lambda)))
g4 = ggplot(df,aes(x=lambda,y=theta,col=variable))+
geom_line()+
theme(legend.position="bottom")+
scale_x_log10()
EN_cv <- cv.glmnet(...) # A COMPLETER
best_lambda <-EN_cv$lambda.min
g4=g4 + geom_vline(xintercept = best_lambda,linetype="dotted",
color = "red")
ggplotly(g4)